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ABSTRACT 

We model the non-linear ideal magnetohydrodynamics of poloidal magnetic fields in neutron stars 
in general relativity assuming a polytropic equation of state. We identify familiar hydromagnetic 
modes, in particular the 'sausage/varicose' mode and 'kink' instability inherent to poloidal magnetic 
fields. The evolution is dominated by the kink instability, which causes a cataclysmic reconfiguration 
of the magnetic field. The system subsequently evolves to new, non-axisymmetric, quasi-equilibrium 
end-states. The existence of this branch of stable quasi-equilibria may have consequences for magnetar 
physics, including flare generation mechanisms and interpretations of quasi-periodic oscillations. 
Subject headings: Magnetohydrodynamics — Stars: magnetars — Stars: magnetic field — Stars: 
neutron 



1. INTRODUCTION 

Neutron stars harbour the strongest known magnetic 
fields in Nature. It is therefore not surprising that the 
magnetic field plays a key role in many aspects of neutron 
star physics. Most of what is known about neutron star 
magnetism comes from the inferred dipole moment as- 
sociated with the electromagnetic stellar spin-down. In 
systems with prominent thermal emission, the surface 
magnetic field strength can be independently measured 
by its influence on the emission pattern. Remarkably, 
various astronomical observations suggest a systematic 
difference in field strengths b etween different classes of 
neutron stars (see lKaspill2010l for a recent review). Aged 
systems, such as millisecond pulsars and accreting neu- 
tron stars in low-mass X-ray binaries, have relatively 
weak magnetic fields, B ^ 10^ G. Typical radio pulsars 
have B ~ 10^^ G. The highest field strengths are associ- 
ated with soft-gamma-repeaters (SGRs) and anomalous 
X-ray pulsars (AXPs) with magnetic fields B ~ 10^^ G. 

This last class of neutron stars, collectively known 
as magnetars, are the focus of this paper. These 
slowly-spinning objects exhibit high energy emission, 
which is occasionaly punctuated by bursts. In some 
rare occasions, powerful flares have been observed in 
SGRs. According to t he prevailing magnetar model, 
as first put forward b y 'Dunca n fc Thompson I (|1992f ): 
[Thompson fc Duncan I (1995), the energy reservoir 
powering magnetar activity is the energy of the 
magnetic field itself. The unusually high magnetar 
surface temperatures are likely the result of heating 
produced by the Ohmic dec ay of the magnetic field in 
the neutron star crust (e.g. iPons et al. I l2007t ). More 
spectacularly, giant SGR flares could be the mani- 
festation of global magnetic fle ld instabilities and of 
subse quent field reconfiguration ([Thompson fc Duncan I 
11995( 1. Furthermore, there is strong evidence that 
quasi-periodic oscillations (QPOs) observed in the 
light curves of these giant flares are of magneto- 
elastic nature (i.e. the restoring force is provided 
by magnetic and crustal elastic stresses) (e.g. iLevin I 

|lasky®tat.physik.uni-tuebingen.de| 
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netar magnetohydrodynamics is hindered by the lack of 
information about the properties of the interior magnetic 
fleld. Such an understanding is a vital element if we are 
to carry out, for example, magnetar 'asteroseismology' 
using the observed QPOs. 

The aforementioned reasons imply it is of great theo- 
retical and astrophysical interest to model the dynamics 
of magnetic fields in the interior of magnetars (and in 
general of neutron stars). More specifically, three key 
questions have attracted the most attention to date: (i) 
What is the nature of long-term magnetic field equilibria 
in neutron stars? (ii) Do hydromagnetic instabilies play 
a role in phenomena like the bursts and flares seen in 
AXPs and SGRs? (iii) Which physical mechanism(s) de- 
termines the long-term evolution and decay of magnetic 
flelds in neutron stars? 

In this letter we make contact with these flrst 
two questions. We study short-term (i.e. dynami- 
cal timescales) neutron star magnetohydrodynamics 
(MHD), speciflcally addressing the issue of global MHD 
instabilities in neutron stars and resulting magnetic 
fleld equilibria. The last few years have seen a flare of 
activity on this subject, mainly motivated by magnetar 
astrophysics. A large body of work has been de- 
voted to semi-a nalyti cal studies of axisymmetric MHD 
equilibria (e.g. Iloka I E OOI; Yoshida fc Eriguchi 2006; 
lYoshida. Yoshida & Eriguchi 2006; HaskeU et al. 2008; 
Ciolfl aL ,2009.; Ciolfi, Ferrari & Gualticri 201Q,), 
build ing on much earlier work regarding equflibria 
(e.g. iMonaghan 1 119651; [Roxburgh 1 11966! ; iParker 1 1196' " 
and stability analys es (Kruskal fc Schw ajzschild I 195 
Tavlerl llW . IgTS": 'Wrigh t I 119731: iMarkev fc Tavler I 
19731 11974: Flowers & Rud erman 1119771 ). It is only 
recently that global, fully non-linear, MHD nu- 
merical calculations have appeared studying mag- 
netic field stability and equil ib ria in stellar models 
(iGeppert fc Rhein hardTI 120061; iBraithwaite fc Spruit 
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iKiuchi. Yoshida fc Shibata~ll201l[ ). The work presented 
herein fahs into this last category. We study the 
dynamical evolution of purely poloidal magnetic fields in 
neutron stars in general relativity, investigating intrinsic 
instabilities and subsequent equilibria, uncovering for 
the first time non-axisymmetric magnetic field equilibria 
in barotropic stars. 

2. NUMERICAL MODEL 

We study magnetized neutron stars through time evo- 
lution of the ideal MHD equations in general relativ- 
ity under the Cowling approximation. To this end, 
we utilise the three-dimensional GRMHD code THOR 
(|Zink. Schnetter fc Tigliol [20081: iKorobkin efdl [20Tol) 
and her sister GPU code horizon ('Zink I l2011[ ). 

The MHD portion of the code utilizes t he co nserva- 
tive formalism outlined in lAnton et al. I (|2006l ). Hy- 
perbolic divergence cle aning is employed according to 
[Anderson et al. I ()2006D . We treat the exterior of the 
star as a low-density artificial atmosphere, allowing full 
evolution of the magnetic field in this region. The outer 
boundary of our domain, located at approximately 1.4 
times the stellar radius, uses Dirichlet boundary condi- 
tions for the magnetic field. We have tested our evo- 
lutions against this condition and found no discernible 
difference with the nature of the instability or the end- 
state of the magnetic field configuration. 

We employ a Cartesian grid with 120'^ grid points at 
single precision. We have completed numerous resolution 
studies including 90^, 150^ and 200^ grid points, finding 
consistent phenomenology with converging growth-times 
for the instability. Moreover, we have tested the use of 
double precision on these simulations and again found 
consistent phenomenology and timescales. 

2.1. Initial Conditions 

We employ the spectral code lorene[3, which pro- 
duces self-consistent solutio ns of the Einstein-M axwell 
field equations in ideal MHD (jBocquet et al. Ill995[ ). The 
spectral grid is mapped to our Cartesian grid and the sys- 
tem is allowed to evolve. No additional perturbations are 
imposed, relying instead on intrinsic numerical noise. 

2.2. Mode Decomposition 

The modal structure of the instability is extracted by 
performing a Fourier decomposition of various physical 
quantities on a ring in the equ atorial plane. W e compute 
complex weighted averages (|Zink et al. ii2007f) 

Cm{f) r f 0, ^ = 0) e^"^d(/., (1) 

where m — + y"^ — const, lies in the initial equa- 
torial plane of the magnetic field. In equation ([ij, / is 
some quantity in the star chosen to best represent the 
instability. Throughout much of this article we use the 
(j) and z components of the magnetic field, although note 
the velocity and density fields also exhibit the instability, 
however the signal is generally not as clean. 

2.3. Fiducial Neutron Star Model 
^ http://www.lorene.obspm.fr/ 




time [ms] 

Fig. I. — Evolution of C'm{B^) as a function of time. Our fiducial 
model has an average magnetic field of B = 1.3 X 10^^ G implying 
an Alfven crossing time of 5.0 ms. The arrows represent the times 
of the three-dimensional snapshots plotted in figures [2] 

Our fiducial model is a non-rotating star with poly- 
tropic equation of state relating pressure, P, and rest- 
mass density, p, through P = Kp^, with F = 2 and 
K = 100. This gives a stellar model with gravitational 
mass M = 1.3 Mq and equatorial radius R — 12.6 km. 
Our fiducia]_ model has an average field strength inside 
the star oi B — 1.3 x 10^^ G, yielding an Alfven crossing 
time of tA — 5.0 ms. Such models have polar magnetic 
field strengths of Bguri = 8.8 x 10^^ G, which is a factor of 
a few larger than observations of the dipole component 
of the strongest magnetic fields. This is predominantly 
done to reduce computational costs. 

3. RESULTS 

In figure [1] we plot the evolution of to = 1 , . . . , 4 
modes for C„i (B^) for our fiducial model, measured at 
vj — 0.6Wi,, where Wi, is the stellar radius. An instabil- 
ity is present in all modes after the first couple of Alfven 
crossing times. Each mode grows exponentially by ap- 
proximately six orders of magnitude over many subse- 
quent Alfven crossing times. Saturation of the modes 
occurs after about 75 ms, at which point the simulation 
evolves to a pseudo-cquillibrium state. We have evolved 
such simulations to 400 ms — 80 1^, with little variation 
in the equilibrium configuration following the first hun- 
dred or so milliseconds. We discuss the equilibrium state 
in more detail below. 

In the first few milliseconds of the simulation we see 
the TO = 4 mode in figure [1] grow by an order of mag- 
nitude over the other modes. This feature is somewhat 
aligned with our Cartesian grid, which we attribute to a 
finite differencing effect associated with the surface of the 
star. To understand this further we performed resolution 
studies with 150'^ and 200"^ grid points. In such simula- 
tions we see the initial transient that lifts the to = 4 
mode above other modes slowly converging away with 
increasing resolution. 

In figures [2] we present three-dimensional plots of the 
magnetic field at various instances throughout the evo- 
lutior0. The instability acts near the neutral line, which 
is the line where B = around the equatorial plane, lo- 
cated in our simulations at approximately two-thirds of 
the stellar radius. For clarity, we have plotted red field 

^ A movie of the fiducial simulation lasting 400 ms can be viewed 
at http: //www. tat .physik .uni-tuebingeii. de/~ tat/grmhd/ 
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lines seeded near the neutral line in the equatorial plane. 
Also plotted are black field lines seeded on the equato- 
rial plane interior to the neutral line. The blue volume 
rendering is an isopycnic surface of p = 0.37pc, where pc 
is the central rest-mass density, which lies at a radius of 
approximately 50% of the star; well inside the neutral 
line. 

Figure [2^ shows the initial data imported from the 
LORENE spectral code. The domain of our grid is larger 
than that plotted here and field lines are truncated at 
the surface of the star for clarity. Such a configuration is 
known to be unstable to the kink ins tability from local 
linear studies in Newtonian physics ()Markev fc Tavlerl 
Il973f ). Such studies have been confirmed with global , 
linearized, numerical evolutions fLander fc Jones 1 120 111) 
and in full non- linear studies (Braithwaitc 2007). This 
is the first time such a configuration has been studied in 
general relativity. 

Figure shows the evolution after 25 ms = 

5tA- The onset of th e 'sausage' or 'varicose' mode 
(|Markev &: Tavler"1ll973l ). involving changes in the cross- 
sectional area of a flux tube around the neutral line, is 
clearly visible. This is strongest in the m = 4 mode, 
which is a result of the transient excitation at the be- 
ginning of our simulation. Whilst this transient reduces 
with increasing grid resolution, the presence of the vari- 
cose mode is an inherent characteristic of the system (as 
evidenced from its almost constant quasi-equilibra state 
that is resolution independent). 

The varicose mode is visually present in our simula- 
tion for approximately eight or nine Alfven crossing times 
before the 'kink' instability appears and begins to domi- 
nate. In figure [2]: we show the onset of this mode acting 
perpendicular to the gravitational field, in accordance 
with the prediction of Markey & Tayler (1973). This fig- 
ure is after 50 ms — 10 1 a, at which point we can still see 
the presence of the varicose mode; i.e. the red field lines 
still exhibit varying cross-sectional area around the star. 
It is worth noting that this mode will have been excited 
from the beginning of the simulation (as seen in figure[T]), 
however it is only now that the exponential growth has 
reached a point at which it is visually obvious. In many 
ways this represents the non-linear development of the in- 
stability where the change in field structure is of similar 
order to the background field. It is also worth mention- 
ing here that our choice of the Cowling approximation 
is well justified considering the dominant instability acts 
on equipotential surfaces. 

The modal analysis presented in figure [T] does not dis- 
tinguish between the varicose or kink modes, implying 
exponential growth of these quantities does not necessar- 
ily indicate an instability in one or the other mode. In- 
deed, figuresOD and[5}: suggest that, if the varicose mode 
is unstable it grows slower than the kink mode. It is the 
kink mode that causes a cataclysmic reconfiguration of 
the topology of the magnetic field. The linear analysis re- 
garding the st ability of the varicose mode is also unclear. 
iTavler I (|1957[ ) showed the varicose mode is unstable in 
cylindrical geometry, however utilising a system that is 
not eq uivalent to a full neutron star. iLander fc Jones I 
(|2011f ) looked at the stability of poloidal fields in neu- 
tron stars by way of numerical evolutions of the lin- 
earized equations, however did not identify the differ- 
ence between the varicose and kink modes. Summarily, 



we believe the stability of the varicose mode to be an 
open question. 

Figure [2ji shows the simulation after 195 ms = 39 Ia- 
This is a typical snapshot many Alfven timescales after 
the non-linear saturation of the unstable modes. We dis- 
cuss the end-state in more detail in section [221 




Fig. 2. — Time evolution of fiducial model with average magnetic 
field B = 12.8 X lO'^^ G, corresponding to an Alfven wave crossing 
time of 5ms. The figures are a) t = ms, b) t = 25 ms, c) t = 50 ms 
and d) t = 195 ms. To more clearly visualise the instability, the red 
field lines are seeded on the equatorial plane close to the neutral 
line, and the black field lines are seeded on the equatorial plane 
interior to the neutral line. The volume rendering is an isopycnic 
surface at 37% of the central rest-mass density, shown to provide 
contrast with the field lines. 



3.1. Instability Growth- Times 

We are interested in extracting instability growth- 
times for the individual modes. We define the growth- 
time for each mode as 



Aln[C™(/)]' 

which we measure during the near exponential growth of 
the instability. We have performed multiple simulations 
with different magnetic field strengths to examine the 
scaling behaviour of the instability. In figure [3] we plot 
the growth-times as a function of mode number for each 
model. The scaling behaviour is consistent with what 
is anticipated, i.e. models with stronger magnetic fields 
and hence faster Alfven velocities have faster growing 
instabilities. This scales almost linearly;_the model with 
average field strength B15 = 4.4 (where B15 = B/IO^'^ G) 
has an Alfven crossing time of tA — 13.9 ms and growth- 
time of approximately t 12 ms, whereas the crossing 
time for the B15 = 25.6 model is tA — 2.4 ms and growth- 
time r ^ 2 ms. 

Whilst this coarse measurement of growth-times scales 
correctly with magnetic field strength, for some of our 
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Fig. 3. — Growth-time, t, as a function of mode number, m, 
for models with average magnetic field strengths ranging from 
Bis = 4.4 - 25.6, where B15 = G. We plot error bars on 

the fiducial, S15 = 12.8 simulation, which are derived from deter- 
mining the growth-time scale, equation l(2]l, over different regions 
of the instability. 

models the individual modes do not scale in the manner 
nece ssarily expecte d. iKruskal fc Schwarzschild I (|1954[ ) 
and iTavler I ()1957[ ) showed the instability grows faster 
for increasing m (our m corresponds to their wavenum- 
ber k) , akin to what is seen for our two simulations with 
the weakest magnetic fields, i.e. models with B15 = 4.4 
and 8.5 in figure |3l This fur ther agrees with the analy- 
sis of iLander fc Jonesi (j2011[ ) , who showed growth-times 
converging soon after the m = 4 mode. We do note that 
our simulations with stronger magnetic fields show dif- 
ferent behaviour in that the growth-time as a function of 
mode number is almost constant. 

We provide three possibilities for the cause of the afore- 
mentioned discrepancy. 1) We are possibly seeing here 
mode coupling between the lower order modes, which 
may scale with magnetic field strength. A more detailed 
analysis of this is required. 2) We are extracting a su- 
perposition of both the varicose and kink modes, which 
potentially have significantly different growth-times for 
the same mode number. 3) Our mode signal is 'dirty', 
in that when evaluating equation ([2|), we pick a region of 
the exponential growth and perform a linear regression 
(in log space) to determine the growth-time. Choosing a 
different region of the exponential growth can give signif- 
icantly varied growth-times. This is displayed by way of 
error bars on the fiducial model in figure |31 which are of 
similar order to the size of variation one expects between 
the modes. We expect this to be an inherent difficulty 
with extracting such growth timescales from non-linear 
simulations. 

3.2. Equilibrium Configurations 

We now return to figure [2ji, particularly concerning 
our equilibrium configurations. In some ways, figure [2Ji 
is consistent with the 'twisted- torus' configurations seen 
in the non- linear evolutions of iBraithwaite fc Nordlund I 
(|2006') and IBraithwaite * (2009), and in the semi-analytic 
equilibriur n derivations of Yoshida. Yoshida & Eriguchi 
([20061) and lCiolfi al. I (|2009t ). This figure shows about 
half of the star is well approximated by a twisted-torus, 
however the remainder of the star also exhibits large non- 
axisymmetric structures. This is typical of the evolu- 
tion following the saturation of the instability. Whilst 



IBraithwaite I ()2008l ) has also found non-axisymmetric 
equilibria as a result of various non-linear simulations, 
we believe ours are the first such simulations to do so 
with a barotropic equation of state. This therefore ex- 
tends existing results in the literature (e.g. iCiolfi et al. I 
l2009f ). in that we provide evidence for a new branch of 
stable, equilibrium solutions with a barotropic equation 
of state. 

A way of classifying variation from the initial state is 
through calculating the relative energy in the toroidal 
and poloidal components of the magnetic field. The ini- 
tial field obviously contains Ep/E — 1.0, where Ep and E 
are respectively the poloidal and total magnetic field en- 
ergy's. Throughout the evolution of the instability we see 
a transfer of energy between the toroidal and poloidal en- 
ergies, although it is interesting to note that this does not 
take place until we visually see the onset of the kink in- 
stability (i.e. '^SOms). Our equilibrium configuration is 
then characterised by a magnetic field with Ep/E ~ 0.75, 
which is similar to the upper-end of the stability window 
found in stably stratified models of IBraithwaite I ([2009). 

It can further be seen from figure [2ji and the full evolu- 
tion of the system that, despite large-scale restructuring 
of the interior field, the center of the star is still threaded 
by a dominantly poloidal component. Moreover, the 
magnetic field exterior to the star is also dominantly 
poloidal, with little visual change to its initial structure. 
This is also important for stabilit y issues associated with 
poloidal fields, particularly the iFlowers &: Ruderman I 
([1977., ') instability which relies on minimising the energy 
in the exterior region of the star (see recent result of 
iMarchant, Reissenegger fc Akgiin II2O1O0 . 

4. CONCLUSION 

We have presented the first, non-linear GRMHD simu- 
lations of hydromagnetic instabilities inherent to purely 
poloidal magnetic fields. In particular, we have evolved 
initially self-consistent solutions of the Einstein-Maxwell 
field equations under the assumption of ideal MHD and 
the Cowling approximation. Purely poloidal magnetic 
field configurations are particularly vulnerable to two 
modes ~ the varicose/sausage and kink modes. The vari- 
cose mode, a cross-sectional change in a toroidal tube of 
magnetic field around the neutral line, was initially ex- 
cited in our simulations due to the Cartesian grid. This 
mode evolved to a pseudo-equilibrium (figure [2}d) before 
it was swamped by the kink instability (figure [2j;), which 
acts perpendicular to the gravitational field. The kink 
instability saturated after approximately 15 Alfven cross- 
ing times, at which point the system evolved towards a 
new, non-axisymmetric equilibrium (figure [2}i). 

The equilibrium end-states achieved from our 
simulations partially resemble twisted-torus con- 
figurations seen from o ther n on-linear ey olutions 
(IBraithwaite fc Nordlundl 120061: IBraithwaite! [2008). 
The center of the star is threaded by a dominantly 
poloidal component, which is consistent with the in- 
stability dominating near the neutral line of the star. 
In some regions of the star we see the presence of 
twisted fiux tubes, with other regions displaying highly 
non-axisymmetric portions of the magnetic field. These 
are the first non-axisymmetric equilibria attained with 
a barotropic equation of state. An interesting extension 
to this work is to study the evolution of the exterior 
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magnetic field as a result of the interior rearrangement. 
Observations of the spin down of neutron stars only 
infers the dipole component of the field, implying the 
presence of extra degrees of freedom can only be inferred 
at this stage through detailed modeling. However, a 
better understanding of the dynamics of the stellar 
surface in our simulations is required to accomplish this 
goal. 

There are many further extensions of this work 
worth pursuing towards more realistic neutron star 
models. Such models should feature composi- 
tion stratification, which is likely to have an 
impact on MHD equilibrium configurations (e.g. 
iReisenegger I [200910 . Further improvements can in- 
clude crustal & magnetospheric physics and study- 
ing the effects of superfiuidity/superconduct ivity (e.g. 
iGlampedakis. Andersson fc Samuelsson II201H ). 

The end-state quasi-equilibria models discov- 
ered herein add to the possible solution space of 
barotropic equilibria. Axisymmetric equilibrium 



models have been explor e d by iCiolfi et al. I (|2009f) : 
iCiolfi. Ferrari fc Gualtieril (|2010D . which would be in- 
teresting to test for stability. Indeed, any pseudo-stable 
equilibrium could serve as possible starting points 
for studying the secular magnetic field evolution in 
magnetars and its connection to activity seen in these 
objects. Moreover, the theoretical modeling of magnetar 
QPOs (which has so far been done using symmetric 
magnetic field configurations), should be extended to 
account for the non-axisymmetric configurations akin to 
those encountered here. We hope to address some of 
these exciting issues in the near future. 
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are grateful to Jerome Novak and Erik Schnetter for as- 
sistance with Lorene, and particularly grateful to Jon 
Braithwaite for useful discussions. 
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